Numerical test of the theory of pseudo-diffusive transmission at the Dirac point of a 

photonic band structure 
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It has recently been predicted that a conical singularity (= Dirac point) in the band structure 
of a photonic crystal produces an unusual 1/L scaling of the photon flux transmitted through a 
slab of thickness L. This inverse-linear scaling is unusual, because it is characteristic of radiative 
transport via diffusion modes through a disordered medium — while here it appears for propagation 
of Bloch modes in an ideal crystal without any disorder. We present a quantitative numerical test 
of the predicted scaling, by calculating the scattering of transverse-electric (TE) modes by a two- 
dimensional triangular lattice of dielectric rods in air. We verify the 1/L scaling and show that the 
slope differs by less than 10% from the value predicted for maximal coupling of the Bloch modes in 
the photonic crystal to the plane waves in free space. 



PACS numbers: 42.25.Bs, 42.25.Gy, 42.70.Qs 
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Two-dimensional photonic crystals with a triangular 
lattice (such as shown in Fig. have been studied ex- 
tensively P, S i i, i i 0, 1 i [O, [S [11, in 
particular because they have a well-developed band gap. 
For frequencies inside this gap the transmission through 
the crystal decays exponentially with the thickness L. 
The band structure has another interesting feature which 
has received much less attention, namely the conical sin- 
gularity that appears at the corner (= iC-point) of the 
hexagonal first Brillouin zone [l| . As indicated in Fig. [21 
at a given wave vector near the if-point two Bloch modes 
are nearly degenerate in frequency. The envelopes of the 
Bloch modes satisfy a pair of coupled differential equa- 
tions that have the same form as the Dirac equation of 
relativistic quantum mechanics p^ . Hence the name 
"Dirac point" given to the conical singularity. The es- 
sential difference between a band gap and a Dirac point 
is that the density of states is zero for a finite frequency 
interval in the former case, but only at a single frequency 
in the latter case. 

Motivated by an electronic analogue (graphene p^). 
Bazaliy and the authors [lB| have recently predicted a 



FIG. 1: Top view of a two-dimensional photonic crystal 
formed by dielectric rods on a triangular lattice in the x — y 
plane, aligned along the z-direction. The lattice constant a 
(centre-to-centre distance of the rods) is indicated. We cal- 
culate the transmission through the slab of thickness L of 
radiation incident near the A"-point of the photonic crystal, 
and find that it scales as 1/L. 
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FIG. 2: Electromagnetic band structure of the photonic crys- 
tal shown in Fig. [TJ calculated for a dielectric constant e — 14 
inside the rods and e = 1 (air, speed of light c) outside the 
rods. (We used the mpb software package for this type of cal- 
culation [isl].) The rods (radius r — 0.27 a) occupy a fraction 
/ — 0.26 of space in the crystal. The bands are shown for the 
case that the magnetic field is parallel to the rods (TE modes). 
The arrow points to the conical singularity (Dirac point) and 
the dashed line shows the dispersion relation in free space. 
The first Brillouin zone is drawn in the inset. (Note that the 
T — Ad direction is perpendicular to the a; = interface of the 
photonic crystal, for the orientation of Fig. [T]) 



new signature of the conical singularity: near the Dirac 
point the photon flux I transmitted through a slab of 
photonic crystal is predicted to scale as 1/L with the 
thickness L of the slab. The 1/L scaling is called 
"pseudo-diffusive" due to its reminiscence of diffusion 
through a disordered medium — although here it appears 
for Bloch modes in the absence of any disorder inside the 
photonic crystal. 

More quantitatively, the prediction of Ref. [lB| is that 
at the Dirac point 



LoFo-^, 0<ro<l/7r, 
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FIG. 3: Equifrequency contours for the photonic crystal of 
Fig- m calculated for ui = 2.89 c/a = 0.95uid- The contours 
(thick solid lines) are centered at the corners of the first Bril- 
louin zone, and are approximately circular with a slight trigo- 
nal distortion. The dashed circle is the equifrequency contour 
in free space, at the same u. An incident plane wave at an an- 
gle 6 (dashed arrow) is coupled to Bloch modes in the crystal 
with the same wave vector component ky (solid arrow shows 
wave vector of the envelope field). When uj —* ljd, the radius 
of the equifrequency contours shrinks to zero and the inci- 
dent plane wave can only couple to evanescent (exponentially 
decaying) Bloch modes. 



with /q the incident photon current per transverse mode 
and / the transmitted photon flux (= transmitted photon 
current per unit width). The coefficient Fq that deter- 
mines the slope of the 1/L scaling depends on the cou- 
pHng strength of the Bloch modes inside the photonic 
crystal to the plane waves outside. For maximal coupling 
one has Fq = I/tt 0,(3. 

It is the purpose of this paper to test the prediction 
of Rcf. 16] quantitatively, by means of a numerical so- 
lution of the scattering problem. (An independent test 
in Ref. [13] provides only a qualitative comparison.) By 
means of an exact solution of Maxwell's equations we 
can test how well the Dirac equation used in Ref. [l^ de- 
scribes the scattering near the Dirac point. Furthermore, 
we can determine the slope Fg — which is beyond the 
reach of the Dirac equation and was left undetermined in 
Ref. [3. 

We solve the scattering problem in the geometry of 
Fig.[T]for the parameters listed in Fig.O The transmitted 
photon flux for a given incident plane wave oc g^^^x+ikyy 

is calculated as a function of frequency w = c^J fc^ -|- ky 

for a given thickness L of the crystal. (The transverse 



width is infinite in the calculation.) We use the finite- 
difference time-domain method [Toj . as implemented in 
the MEEP software package |20l. 

To make contact with Ref. p6| we first extract from 
Fig. m the parameters lud — 3.05 c/a, vd — 0.369 c that 
characterise the conical singularity in the band structure, 

Slu = Lu — LU]j = VD\Sk\. (2) 

Here Sk — k — K is the displacement of the wave vector k 



from the X-point, with wave vector K = |7ra~^(\/3, 1). 
The velocity w_d is the group velocity of Bloch modes at 
frequencies near the frequency lod of the Dirac point. A 
given Sky corresponds to an angle of incidence 



= arcsm 



-(Ky+Sky) 



(3) 



In particular, Sky = and lu = ujd correspond to = 
arcsin (27rc/3wz)a) = 6*0. For our parameters 9o — 43°. 

As indicated in Fig. [31 an incident plane wave couples 
to Bloch modes in the photonic crystal with the same 
ky. Propagating envelope modes have wave vector on 
the equifrequency contour centered at a X-point. As the 
frequency uj approaches the Dirac frequency ujd, the ra- 
dius of the equifrequency contour shrinks to zero, and the 
incident plane wave can only couple to evanescent modes. 
These decay exponentially away from the interface, with 
a decay length cx l/](5fcj^] which becomes infinitely long 
at the ii'-point. 

The crucial difference between transmission at the 
Dirac frequency and inside a band gap is this: In both 
cases, the photonic crystal supports only evanescent 
Bloch modes, but inside the band gap the decay length 
as a function of angle of incidence has a finite maximum 
value — while at the Dirac frequency the maximum decay 
length is infinite. As a consequence, angular averaging of 
the transmitted intensity over some narrow range of in- 
cident angles around Oq gives an exponentially decaying 
transmission inside the band gap, but only an algebraic 
1/L decay at the Dirac frequency [I6I] . 

For a quantitative description of this scaling behavior 
we need to consider the coupling strength of the Bloch 
modes inside the crystal to the plane waves outside. The 
transfer matrix of the interface at x = and x = L, 
which determines this coupling, is characterised by two 
parameters (3 and 7. These parameters enter into the 
expression for the transmission probability T(6ky, 6lo), 
which is defined as the ratio of transmitted to incident 
photon flux for an incident plane wave [frequency uj = 
lud + Suj and angle of incidence 9 related to Sky by Eq. 
(I2D]. The result is [l^] 
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with a = {5ui/vdY — Sky. 

We extract the two interface parameters 

(3 = -0.094, 7 = -0.133 



(5) 



from the T versus Suj dependence at Sky = 0, plotted in 
Fig. m In the same figure we show that the Sky depen- 
dence of these parameters is weak for Skya ^ 1, as was 




FIG. 4: Transmission probability through the slab of pho- 
tonic crystal of thickness L — 8-\/3a. The data points are the 
numerical results, the curves are calculated from Eq. Q with 
the interface parameters of Eq. (O . The vertical dashed line 
indicates the Dirac frequency ujd ■ This plot is for a single inci- 
dent plane wave with Sky = (open data points, solid curve) 



and Sky 



-(7r/30)a (filled data points, dotted curve) 
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FIG. 6: Transmitted flux Jmin at the minimum near the Dirac 
point versus the thickness L of the slab. Open data points are 
for A = (7r/30)a~^, filled data points are for A = (7r/15)a~^. 
The solid and dashed lines show the analytical prediction from 
Eq. with the interface parameters of Eq. © . 
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FIG. 7: Minimal transmitted flux versus slab thickness for the 
four data sets tabulated in Table|T] We took A = (7r/15)a~^ 
in each case. The analytical result for maximal coupling (/3 = 
7 = 0) is indicated by the solid curve. 



FIG. 5: Transmitted flux ((6)1 for L = 13\/3a. Data points 
are the numerical results, curves are calculated from Eq. 
The vertical dashed line indicates the Dirac frequency ujd- 
This plot is for a range \5ky\ < A of incident wave vectors, 
with A — (7r/30)a~^ for the open data points and solid curve; 
A = (7r/15)a~^ for the flUed data points and dotted curve. 

assumed in Ref . , since the same set of parameters JS]) 
also describes the Sco dependence of T at nonzero Sky. 

To test for the 1/L scaling we need to consider a range 
—A < Sky < A of incident transverse wave vectors. (This 
corresponds to a range A6 ~ 2cA/a;j) cos On of incident 
angles centered at Oq.) According to Ref. the 1/L 
scaling is reached when i > 1/A. We calculate the trans- 
mitted photon flux I{lu) in this range of wave vectors, 

/M=/o/ -^T{6ky,SLj^uj-ujD). (6) 
As shown in Fig. [51 we find a strong dependence of / 



on the range of wave vectors A away from the Dirac 
frequency — but not at the Dirac frequency, where the 
transmitted flux reaches a minimum Imin which is A in- 
dependent for A > 

In Fig. [S] we plot the L dependence of /min on a 
double-logarithmic scale. For L A^^ the predicted 
1/L scaling of Eq. ([T]) is obtained, with a coefficient 
Fq — 0.30. This coefficient is just 6% smaller than the 
value Fq = I/tt reached for maximal coupling of Bloch 
modes and plane waves at the interfaces between the pho- 
tonic crystal and free space. 

To investigate how generic these results are, we have re- 



The frequency WYnin of the transmission minimum is slightly 
offset from the Dirac frequency uijj , but the relative offset is small 
and vanishes with increasing L: |ix)i„in — uJd 
We have checked that it makes no difference for the 1/L scaling 
whether we calculate the transmitted flux at uJmin or at wjj. 
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TABLE I: Parameters representing four different triangular 
lattice photonic crystals. Symbols on the left correspond to 
the data points in Fig. [7] 

peated the calculation for different values of the dielectric 
constant e of the rods and for different filling fractions / 
(related to the radius r of the rods by / = 27rr^/V3a^). 
The parameters corresponding to the four sets of data 
are tabulated in Table HI In Fig. [7] we show the L de- 
pendence of the minimal transmitted flux for each data 
set. In each case we find 1/L scaling with a slope Fq that 
remains within 8% of the maximal value Fq = I/tt. 
In conclusion, we have presented a quantitative numer- 



ical test of the applicability of the Dirac equation [lj| to 
a photonic crystal with a conical singularity in the band 
structure. The numerical results are in good agreement 
with the analytical predictions flG*! for the transmission 
through a finite slab. In particular, our numerical cal- 
culation demonstrates the 1/L scaling of the transmitted 
photon flux with a slope that is close to the value for max- 
imal coupling at the interface with free space. This find- 
ing implies that transmission experiments can be used to 
search for intrinsic properties of the Dirac point in the 
band structure, not hindered by a weak coupling to the 
outside. 



Acknowledgments 

We have benefited from discussions with M. de Dood. 
This research was supported by the Dutch Science Foun- 
dation NWO/FOM. 



[1] M. Plihal and A. A. Maradudin, Phys. Rev. B 44, 8565 

(1991) . 

[2] P. R. Villeneuve and M. Piche, Phys. Rev. B 46, 4969 

(1992) . 

[3] K. Sakoda, Phys. Rev. B 52, 8992 (1995). 

[4] N. Susa, J. Appl. Phys. 91, 3501 (1995). 

[5] J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Pho- 
tonic Crystals: Molding the Flow of Light (Princeton 
University Press, Princeton, NJ, 1995). 

[6] D. Cassagne, C. Jouanin, and D. Bertho, Phys. Rev. B 
53, 7134 (1996). 

[7] M. Notomi, Phys. Rev. B 62, 10696 (2000). 

[8] S. Foteinopoulou and C. M. Soukoulis, Phys. Rev. B 67, 
235107, (2003). 

[9] A. Berrier, M. Mulot, M. Swillo, M. Qiu, L. Thylen, 
A. Talneau, and S. Anand, Phys. Rev. Lett. 93, 073902 
(2004). 

[10] K. Guven, K. Aydin, K. B. Alici, C. M. Soukoulis, and 
E. Ozbay, Phys. Rev. B 70, 205125 (2004). 

[11] P. V. Parimi, W. T. Lu, P. Vodo, J. Sokoloff, J. S. Derov, 
and S. Sridhar, Phys. Rev. Lett. 92, 127401 (2004). 

[12] R. Moussa, S. Foteinopoulou, Lei Zhang, G. Tuttle, K. 



Guven, E. Ozbay, and C. M. Soukoulis, Phys. Rev. B 71, 
085106 (2005). 

[13] R. Gajic, R. Meisels, F. Kuchar, and K. Hingerl, Phys. 

Rev. B 73, 165310 (2006). 
[14] F. D. M. Haldane and S. Raghu, Phys. Rev. Lett. 

100, 013904 (2008); S. Raghu and F. D. M. Haldane, 

cond-mat/0602501 
[15] J. Tworzydlo, B. Trauzettel, M. Titov, A. Rycerz, and 

C. W. J. Beenakker, Phys. Rev. Lett. 96, 246802 (2006). 
[16] R. A. Sepkhanov, Ya. B. Bazaliy, and C. W. J. 

Beenakker, Phys. Rev. A 75, 063813 (2007). 
[17] X. Zhang, Phys. Lett. A 372, 3512 (2008). 
[18] S. G. Johnson and J. D. Joannopoulos, Optics Express 

8, 173 (2001). 

[19] A. Taflove and S. C. Hagness, Computational Electro- 
dynamics: The Finite-Difference Time-Domain Method 
(Artech House, 2005). 

[20] A. Farjadpour, D. Roundy, A. Rodriguez, M. Ibanescu, 
P. Bermel, J. D. Joannopoulos, S. G. Johnson, and G. 
W. Burr, Optics Letters 31, 2972 (2006). 



